library(leaflet)
library(tidyverse)
library(sf)
library(cancensus)
library(sp)
library(ggplot2)Visualizations of OC Transpo Data Against Canada Census Data
About the Data
The data used in this analysis is sourced from:
- OC Transpo’s General Transit Feed Specification (GTFS), which provides detailed information on Ottawa’s public transportation system, including bus stop locations, route frequencies, and trip schedules1.
- Statistics Canada’s Census Data, accessed via the
cancensusR package, which provides demographic and socioeconomic information at the dissemination area level, including population density, income levels, and commuting patterns.
Libraries
Data Acquisition
The following code loads and organizes the original data sets. The vectors specified in the get_census() function refer to columns grabbed from cancensus. In this case, v_CA21_1004 refers to average total household income, and v_CA21_7635 refers to number of respondents who reported using a car, truck, or van as main mode of transportation to their workplace. Note that dataset="CA21" specifies the data comes from the 2021 Canadian Census.
# Load file of all OC Transpo bus stops
stops <- read_csv("stops.txt")
# Load 2021 census data (for Ottawa) with average income and car as main mode of transportation
ottawa_census_sf <- get_census(dataset = "CA21",
regions = list(CSD = "3506008"),
vectors = c("v_CA21_1004", "v_CA21_7635"),
level = "DA",
use_cache = TRUE, geo_format = 'sf')
# Transform coordinate system (CRS) to latitude/longitude
ottawa_census_sf$geometry <- st_transform(ottawa_census_sf$geometry, CRS("+init=epsg:4326"))Visualizing Ottawa’s Bus Stops and Census Dissemination Areas
leaflet() |>
addTiles() |>
addPolygons(data = ottawa_census_sf$geometry, color = "violet", popup = ottawa_census_sf$GeoUID) |>
addCircleMarkers(data = stops,
lng = ~stop_lon,
lat = ~stop_lat,
radius = 4,
popup = ~stop_name,
fillOpacity = 0.7,
color = 'cornflowerblue')
This interactive map displays all census dissemination areas (DAs) in Ottawa, outlined in "violet", and OC Transpo bus stops, marked in "cornflower blue". Each DA is labeled with its unique identifier (GeoUID), while each bus stop is labeled with its stop name.
This is a basic map, but there are many ways to enhance it with more meaningful insights. For example, I am interested in analyzing the number of bus stops per census area, the average income per polygon, and other transit-related patterns. Below is the code that incorporates these elements into the map.
Mapping Transit Accessibility and Socioeconomic Factors in Ottawa
# Transform stops coordinate system to the same one as ottawa_census_sf to facilitate merging
stops_sf <- st_as_sf(stops, coords = c("stop_lon", "stop_lat"), crs=CRS("+init=epsg:4326"))
stops_sf <- st_transform(stops_sf$geometry, CRS("+init=epsg:4326"))
# Count total bus stops in each census polygon
stops_per_polygon <- st_intersects(ottawa_census_sf, stops_sf) |>
lengths()
# Add as a new column to census data
ottawa_census_sf$stop_count <- stops_per_polygon
# Add bus stop ratio column
ottawa_census_sf$stop_ratio <- ottawa_census_sf$stop_count / ottawa_census_sf$`Shape Area`
# Add average income column
ottawa_census_sf$avg_income <- as.numeric(ottawa_census_sf$`v_CA21_1004: Average total income in 2020 ($)`)
# Add vehicle use ratio per polygon column
ottawa_census_sf$vehicle_ratio <- ottawa_census_sf$`v_CA21_7635: Car, truck or van` / ottawa_census_sf$Population
# The map with avg income, bus stop count, and bus stop layer
leaflet() |>
addTiles() |>
# Layer for bus stop count
addPolygons(data = ottawa_census_sf,
fillColor = ~colorNumeric("Blues", stop_count)(stop_count),
color = "black",
fillOpacity = 0.7,
group = "Total Stops",
popup = ~paste("Avg Total Income: $", format(avg_income, big.mark=","),
"<br />", "Total Stops: ", format(stop_count, big.mark=","),
"<br />", "Vehicle Use Ratio (%): ", format(round(vehicle_ratio,3)*100, big.mark=","),
"<br />", "Stop per sq km: ", format(round(stop_ratio,3), big.mark=","),
"<br />", "Area (sq km): ", format(`Shape Area`, big.mark=","))) |>
# Layer for Stop per sq km
addPolygons(data = ottawa_census_sf,
fillColor = ~colorNumeric("Oranges", vehicle_ratio)(vehicle_ratio),
color = "black",
fillOpacity = 0.7,
group = "Vehicle Use Ratio (%)",
popup = ~paste("Avg Total Income: $", format(avg_income, big.mark=","),
"<br />", "Total Stops: ", format(stop_count, big.mark=","),
"<br />", "Vehicle Use Ratio (%): ", format(round(vehicle_ratio,3)*100, big.mark=","),
"<br />", "Stop per sq km: ", format(round(stop_ratio,3), big.mark=","),
"<br />", "Area (sq km): ", format(`Shape Area`, big.mark=","))) |>
# Layer for Avg Total Income (overlay vehicle ratio and bus stop count as well)
addPolygons(data = ottawa_census_sf,
fillColor = ~colorNumeric("Reds", avg_income)(avg_income),
color = "black",
fillOpacity = 0.7,
group = "Avg Total Income",
popup = ~paste("Avg Total Income: $", format(avg_income, big.mark=","),
"<br />", "Total Stops: ", format(stop_count, big.mark=","),
"<br />", "Vehicle Use Ratio (%): ", format(round(vehicle_ratio,3)*100, big.mark=","),
"<br />", "Stop per sq km: ", format(round(stop_ratio,3), big.mark=","),
"<br />", "Area (sq km): ", format(`Shape Area`, big.mark=","))) |>
addPolygons(data = ottawa_census_sf,
fillColor = ~colorNumeric("Blues", stop_ratio, domain = c(0,max(stop_ratio)))(stop_ratio),
color = "black",
fillOpacity = 0.7,
group = "Stop per sq km",
popup = ~paste("Avg Total Income: $", format(avg_income, big.mark=","),
"<br />", "Total Stops: ", format(stop_count, big.mark=","),
"<br />", "Vehicle Use Ratio (%): ", format(round(vehicle_ratio,3)*100, big.mark=","),
"<br />", "Stop per sq km: ", format(round(stop_ratio,3), big.mark=","),
"<br />", "Area (sq km): ", format(`Shape Area`, big.mark=","))) |>
# Layer for Area (sq km)
addPolygons(data = ottawa_census_sf,
color = "black",
group = "Area (sq km)",
popup = ottawa_census_sf$`Shape Area`) |>
# Layer control to toggle between layers
addLayersControl(
baseGroups = c("Total Stops", "Avg Total Income", "Vehicle Use Ratio (%)", "Stop per sq km"),
options = layersControlOptions(collapsed = FALSE)
)